With domcol module, you can visualize complex functions by domain coloring. You can also draw isolines for the modulus (magnitude) and the phase (argument) of a complex array.
image visualizes a complex function f over a rectangle defined by the lower left point a and the upper right point b. rtpalette draws a color palette and labels.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
image(f, a, b);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(pa, pb);
You can set a scaling factor s for the modulus-axis. The factor must be a positive number and the default value is 1.
The functions and are used by default.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
image(f, a, b, s);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, pa, pb);
You can draw iso-modulus lines with isomodulus and iso-phase lines with isophase. Calculating function values at each point in advance is recommended. rtpalette returns tick values except and . You can use the tick values to draw iso-modulus lines.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
pair[][] z = seq(f, a, b);
int[] H = image(z, a, b, s);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
real[] rs = rtpalette(s, pa, pb);
isomodulus(z, a, b, rs);
isophase(z, a, b);
draw(box(a, b));
You can visualize a domain of a circle of radius rmax with scaling factor c.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
real rmax = infinity, c = 3;
pair[][] z = seq(f, rmax, c=c);
real s = 3;
int[] H = image(z, s);
real[] rs = rtpalette(s, H);
isomodulus(z, rs);
isophase(z);
raxis(rmax, c);
taxis();
You can set a modulus range with the minimum and the maximum value. The default range is (0, infinity).
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
pair rlim = (0.3, 30);
image(f, a, b, s, rlim);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, pa, pb, rlim);
You can also set an approximate modulus range with autoscale=true, or the lower and the upper quantiles. image returns an integer array, which contains the information of the specified range. You can get the range with modulus_range.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
pair rlim = (0.01, 0.9); // trim lower 1% and upper 10%
int[] H = image(f, a, b, s, autoscale=true, rlim);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, pa, pb, modulus_range(s, H));
You can reverse a color palette by setting reverse=true in the argument of image. I prefer a visualization with dark at infinity.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
image(f, a, b, s, reverse=true);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, pa, pb, reverse=true);
The palette with Luv color space is used by default. You can use the palette with HSY color space for higher contrast or the OKLab color space for more uniform color. A color palette can be reversed with the negative vale of b in a palette function.
import domcol;
size(15cm);
pen pal(real l, real t) { return palHSY(l, t, a=0.6, b=-0.4); }
// pen pal(real l, real t) { return palLab(l, t, b=-0.5); }
// pen pal(real l, real t) { return palOKLab(l, t, b=-0.5); }
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 3;
image(f, a, b, s, pal);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, pa, pb, pal);
The modulus-axis is divided by nr (default: 8) and the phase-axis is divided by nt (default: 2). The ticks and labels of the modulus-axis or the phase-axis can be turned off with nr=0 or nt=0. To specify the tick values of modulus, you can pass the values (the array of real) to rtpalette.
import domcol;
real[] rs = {0.5, 1, 2};
unitsize(3cm, -7cm);
picture pic;
unitsize(pic, 3cm, 2.5cm);
erase(pic);
rtpalette(pic);
add(pic, (0, 0));
erase(pic);
rtpalette(pic, reverse=true);
add(pic, (1, 0));
erase(pic);
rtpalette(pic, rlim=(infinity, 0));
add(pic, (2, 0));
erase(pic);
rtpalette(pic, rlim=(infinity, 0), reverse=true);
add(pic, (3, 0));
erase(pic);
rtpalette(pic, nr=4, nt=0);
add(pic, (0, 1));
erase(pic);
rtpalette(pic, rs);
add(pic, (1, 1));
erase(pic);
rtpalette(pic, rs, endlabels=false);
add(pic, (2, 1));
erase(pic);
rtpalette(pic, rlim=(0.5, 2), nr=4, tlim=(0, 2), nt=1);
add(pic, (3, 1));
If you set a scaling factor s to 0, image uses quantiles of the modulus.
image returns an integer array, a cumulative weighted histogram, which is used to get quantiles and the range of complex modulus. rtpalette uses the array to draw labels.
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 0;
int[] H = image(f, a, b, s, reverse=true);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(s, H, pa, pb);
You can reuse the previously calculated cumulative weighted histogram with another image.
import domcol;
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
picture pic;
erase(pic);
unitsize(pic, 0.75cm);
pair a = (-5, -5), b = (5, 5);
real s = 0;
int[] H = image(pic, f, a, b, s, reverse=true);
xaxis(pic, Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(pic, Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(pic, box(a, b));
pair a = (-2, -2), b = (3, 3);
draw(pic, box(a, b));
add(pic, (0cm, 0));
erase(pic);
unitsize(pic, 0.75cm * 2);
image(pic, f, a, b, s, H); // use previous H
xaxis(pic, Bottom, a.x, b.x, RightTicks(N=5));
yaxis(pic, Left, a.y, b.y, LeftTicks(N=5));
draw(pic, box(a, b));
pair pa = (b.x + 0.4, a.y);
pair pb = (b.x + 1.2, b.y);
rtpalette(pic, s, H, pa, pb);
add(pic, (9cm - 0.75cm, -0.75cm));
You can also reuse the previously calculated cumulative weghted histogrm with another isomodulus.
import domcol;
picture pic;
erase(pic);
unitsize(pic, 0.75cm);
pair a = (-5, -5), b = (5, 5);
real s = 0;
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair[][] z = seq(f, a, b);
int[] H = image(pic, z, a, b, s, reverse=true);
xaxis(pic, Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(pic, Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(pic, box(a, b));
label(pic, "$f(z)$", ((a.x + b.x) / 2, b.y), 2N);
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
real[] rs = rtpalette(pic, s, H, pa, pb);
isomodulus(pic, z, a, b, rs);
isophase(pic, z, a, b);
add(pic, (9cm, 0));
erase(pic);
unitsize(pic, 0.75cm / 10);
pair a = (-50, -50), b = (50, 50);
pair f(pair z) { return z; }
pair[][] z = seq(f, a, b);
image(pic, z, a, b, s, H); // use previous H
xaxis(pic, Bottom, a.x, b.x, RightTicks(N=2));
yaxis(pic, Left, a.y, b.y, LeftTicks(N=2));
draw(pic, box(a, b));
label(pic, "$z$", ((a.x + b.x) / 2, b.y), 2N);
isomodulus(pic, z, a, b, rs); // use previous rs
isophase(pic, z, a, b);
add(pic, (0, 0));
The image function defined in domcol.asy returns an integer array. The array is the cumulative histogram of a weighted histogram of with a weight of 8 for each point. The specified lower and upper limits are encoded with a weight of 1 and 2 in the histogram. A reverse flag is encoded with a weight 4 at the last bin.
A cumulative distribution function is obtained from the cumulative weighted histogram, which also provides approximate quantiles.
import domcol;
unitsize(8cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
pair a = (-5, -5), b = (5, 5);
real s = 0;
int[] H = image(pic=null, f, a, b, s);
write(" q l r");
for (int i: sequence(8 + 1)) {
real q = 1 - i / 8;
write(
format("%6.3g", q) +
format("%10.3g", quantile(q, H)) + format("%10.3g", modulus(q, s, H)));
}
pair a = (0, 0), b = (1, 1);
pair[] xy = pairs(sequence(H.length) / H.length, H / H[H.length - 1]);
draw(operator --(... xy), red);
xaxis(Label("$l$", MidPoint), a.x, b.x, RightTicks(N=4));
yaxis(Label("$q$", MidPoint), a.y, b.y, LeftTicks(N=4));
string rlabel(real l) {
return (1 - realEpsilon <= l) ? "$\infty$" : format("$%.3g$", rfn(l));
}
xaxis(
Label("$r$", MidPoint), Top, a.x, b.x, LeftTicks(ticklabel=rlabel, N=4));
draw(box(a, b));
q l r
1 1 5.64e+102
0.875 0.977 42.2
0.75 0.971 33.7
0.625 0.964 26.5
0.5 0.952 20
0.375 0.93 13.3
0.25 0.881 7.39
0.125 0.739 2.84
0 0 0
A lightness function maps the infinite interval to the finite interval .
You can use with set_lfn(2) and with set_lfn(0). The default of a lightness function is
import domcol;
size(15cm);
pair f(pair z) { return (z * z - 1) * (z - (2, 1))**2 / (z * z + (2, 2)); }
set_lfn(2);
// set_lfn(0);
pair a = (-5, -5), b = (5, 5);
image(f, a, b);
xaxis(Bottom, a.x, b.x, RightTicks(N=2, n=5));
yaxis(Left, a.y, b.y, LeftTicks(N=2, n=5));
draw(box(a, b));
pair pa = (b.x + 0.5, a.y);
pair pb = (b.x + 2, b.y);
rtpalette(pa, pb);
The function is related to the projection of a unit sphere from the South Pole to the equatorial plane.